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Abstract 

We present a Polar Coordinate Lattice Boltzmann(PCLB) model for compressible flows. A 
method to evaluate the continuum distribution function from the discrete distribution function 
is indicated. Within the model, the temporal and spatial evolutions are treated with via the 
operator-splitting scheme. The temporal evolution is calculated analytically and the convection 
term is solved via a Modified Warming-Beam (MWB) scheme. Within the MWB scheme a suitable 
switch function is introduced. The new PCLB model works not only for subsonic flows but also 
for supersonic flows. It is validated and verified via the following well-known benchmark tests: 
(i) the rotational slip flow, (ii) the stable shock tube problem, (iii) the Richtmyer-Meshkov(RM) 
instability, (iv) the Kelvin-Helmholtz(KH) instability. As an original application, we studied the 
non-equilibrium characteristics of the system around three kinds of interfaces, the shock wave, 
the rarefaction wave and the material interface, for two specific cases. In one of the two cases, the 
material interface is initially perturbed and consequently the RM instability occurs. It is found that, 
the effects of deviating from thermodynamic equilibrium at the material interface differ significantly 
from those at the mechanical interfaces. The coupling effect of molecular motions in different 
degrees of freedom is much more pronounced at the material interface. The initial perturbation at 
the material interface enhances this coupling effect. The amplitude of deviations from equilibrium 
at the shock wave is much higher than those at the rarefaction wave and material interface. Our LB 
results confirm that the temperature increase first initiates the increase of internal kinetic energy 
in the degree of freedom corresponding to the direction of temperature gradient, then increases 
those in other degrees of freedom. By comparing each component of the high-order moments and 
its value in equilibrium, we can draw qualitative information on the actual distribution function. 
These results deepen our understanding of the mechanical and material interfaces from a more 
fundamental point of view, according to the distribution function in continuum velocity space, and 
provide valuable information to improve the physical modeling at a macroscopic level. 

PACS numbers: 47.11.-j, 47.40.-x, 47.55.-t, 05.20.Dd 
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I. INTRODUCTION 



During recent decades the lattice Boltzmann (LB) modeling and simulation have achieved 
great success in various complex flows [l). However, most of these studies were focused on 
nearly incompressible flow, while with increasing the Mach number, the compressibility of 
flows has to be taken into account. Such high speed compressible flows are ubiquitous in 
aerophysics, astrophysics, explosion physics, medical physics and others. Given the great 
importance of shock waves in many fields of physics and engineering, constructing LB models 
for high speed compressible flows has attracted considerable interest since the early days of 
LB research [lj]. 

In 1992 Alexander et al Q formulated a compressible LB model for flows at high Mach 
number via introducing a flexible sound speed. This model works only for nearly isothermal 
compressible systems. In 1999 Yan et al [3] proposed a LB scheme for compressible Euler 
equations. In the years of 1998 and 2003 Sun and his coworker 

0,3 

presented an adaptive LB 

scheme for the two- and three-dimensional systems, respectively. In this model the particle 
velocities vary with the Mach number and internal energy, so that the particle velocities 
are no longer constrained to fixed values. All of those models belong to the standard LB 
framework. However, due mainly to stability problems, to date, applications of LB methods 
to compressible flows remain scanty. 

Besides the standard LB framework, the other way to formulate LB for high speed flows 
is to use the Finite-Difference (FD) scheme to calculate the temporal and spatial derivatives 

n 

of the distribution function. In 1997 Cao et al [6j proposed to use the FD scheme to improve 
the numerical stability and apply nonuniform grids in the LB method. In the past decade, 



Tsutahara, Watari and Kataoka [7 



lit] proposed several nice FDLB models for the Euler and 



Navier-Stokes equations, where the discretizations in the space and in particle velocity are 
separated. In 2005 Xu Q extended the idea to handle binary fluids. However, similar 
to the case of standard LB models, these FDLB schemes only work for subsonic flows. For 
modeling and simulating high speed compressible flows, especially those with shocks, many 
attempts and considerable progress have been achieved fl4- 33] . 

It should be pointed out that, up to now, most of LB models for compressible fluids 
are based on the Cartesian coordinate system. In many cases the flows show divergent, 
convergent, and/or rotational behaviors, for example, in cylindrical or spherical devices. For 



such flow systems, LB models based on polar coordinates, cylindrical coordinates or spherical 
coordinates are more convenient and are less exposed to numerical errors. There have been 
a number of LB methods based on curvilinear coordinates or for axisymmetric cylindrical 
coordinate system. Early in 1992, Nannelli and Succi presented a general framework 
to extend the lattice Boltzmann equation to arbitrary lattice geometries. In this work, a 
finite- volume formulation of LB equation was given. Then some other finite- volume versions 
of the LB method were proposed for irregular meshes {35-39]. In 1997 He and Doolen [loj] 
extended the LB method to apply to general curvilinear coordinate systems via using an 
interpolation-based strategy. In the following year Mei and Shyy [4l| developed a FDLB 
method in body-fitted curvilinear coordinates with non-uniform grids. Later, Halliday et 



al 



421 ] proposed a Polar Coordinate Lattice Boltzmann(PCLB) method for hydrodynamics. 



In 2005 Premnath and Abraham 



431 ] presented a LB model for axisymmetric multiphase 



flows. In this work source terms were added to a two-dimensional standard LB equation for 
multiphase flows such that the emergent dynamics can be transformed into the axisymmetric 
cylindrical coordinate system. But all those LB methods work only for isothermal and nearly 
incompressible flows. In 2010 Asinari et al j^j] formulated a LB scheme to analyze the 
radiative heat transfer problems in a participation medium, but did not take into account 
the effects of fluid flow. In 2011 Watari |45|] formulated a polar coordinate FDLB scheme to 
investigate the rotational slip flow problems in coaxial cylinders. This work presents valuable 
information on the LB application to the cylindrical system. However, this model works also 
only for subsonic compressible flow systems. In the present work we extend the FDLB model 
based on polar coordinates to compressible flow systems with high Mach number so that it 
can be used to simulate flows with shock waves. 

The rest of the paper is structured as follows. In section II we first review the polar FDLB 
model by Watari, then present our contributions for the polar FDLB model. The latter 
includes the operator-splitting scheme, the analytical solution for the temporal evolution 
and the Modified Warming-Beam(MWB) scheme for the convection behavior. Section III is 
for the validation and verification of the new LB model. In section IV we study the non- 
equilibrium characteristics of the system in two special cases related to shock wave passing 
material interfaces. The method to qualitatively recover the actual distribution function is 
illustrated. Section V concludes the present paper. 
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II. POLAR FDLB MODEL 



A. Brief review of Watari model 

Below is a general description of the two-dimensional FDLB thermal model [45] , which is 
applicable to both rectangular cartesian coordinate system and polar coordinate system. The 
evolution of the distribution function f ki with the Bhatanger-Gross-Krook approximation 
j^l reads, 

^+v«-V/ w = — [/«-.£?] (1) 

where fki (fki) 18 the discrete (equilibrium) distribution function; r is the relaxation time 
determining the speed of approaching equilibrium; v/^ is the discrete velocity which will 
be defined below. The original Discrete- Velocity-Model(DVM) by Watari and Tsutahara is 
composed of Nk + 1 groups of discrete velocities. The fc-th group has the size Vk- Each group 
has Ni components distributed in Ni directions. Mathematically, the DVM can be write as: 



^ ^ ^kiot^a ^kix^x ^kiy^yj (2) 



where Vk% x = Vk cos[27r(i — 1)/Nj\, Vk% y = Vk sin[27r(i — 1)/A^], e x and e y are the unit vectors 
in the two-dimensional rectangular cartesian coordinate system, k = 0,1,2,- • • ,A^, and i = 
1,2,- • • ,A^. In this work we discuss the polar coordinate FDLB model for fixed Nk = 4 and 
flexible A^. The sizes of discrete velocities are chosen as v = 0, Vi = 1, v 2 = 2.92, v% = 2.99, 
= 4.49. The sketches of the DVM for cases of Ni = 8, Ni = 16, Ni = 24 are shown in 
Figffl 



It's easy to prove that this DVM with Ni — 8 has at least up to seventh rank isotropy[2l|. 
The macroscopic quantities are defined as 

ki ki 

P u = E fki V ki = E f kiWki ' ( 4 ) 

ki ki 

P E = J2 \fki( v ki - u) • (v w - u) = \fki{v ki - u) • (v w - u). (5) 



2 KlK ' v ' ^2* 

ki ki 



Here p, u (= u r e r + u$eo = u x e x + u y e y ), P (= pT), E{= T / (7 — 1)) are the hydrodynamic 
density, flow velocity, pressure and internal kinetic energy per unit mass, respectively; T is 
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i = 7 

{a) 



i = 13 



i = 19 

(c) 



FIG. 1: Sketches of DVM with fixed N k = 4 and various values of N { . (a) AT* = 8. (b)7V i =16. 
(c)iVi = 24. 

the temperature and 7(= 2) is the specific-heat ratio. Other velocity moments that the local 
equilibrium distribution function has to satisfy are: 

X>7v^ = p(£I + uu), (6) 
Y fki v kiv k iv k i = p[E(u a e^5^ + e a u /3 e 7 <5 7a + e^u^^) + uuu], (7) 
Yl \fki v ki ■ VfciV fei = p{EI + uu) , (8) 

ki 

Yl \fki v ki ■ v ki v ki v ki = p(2E + • u)(El + uu), (9) 

ki 

where I is the unit tensor, v^v^ and uu are double dyadics, v^v^v^ and uuu are triple 
dyadics. 

The equilibrium distribution function is computed by, 

m = pw - ^ + ^) + ^f £ ( l - ^) + - ^) 

VkieVkiirVkiOUsUirUd VkieVkmVkii)Vkie,U £ U n U^U i 

QE 3 24E 4 J 1 j 

where the weighting coefficients are calculated in the following way, 

1 



[B 4 E 4 + B 3 (v 2 k+1 + v£ +2 + v 2 k+3 )E* 



vl{vl-vl +1 )(vl-v 2 k+2 ){vl-vl + zy 
+B 2 (v 2 k+1 v 2 k+2 + v 2 k+2 v 2 k+3 + v 2 k+3 v 2 k+1 )E 2 + B x v 2 k+1 v 2 k+2 v 2 k+3 E\ (11a) 

F = 1/B - (Fx + F 2 + F 3 + F 4 ), (lib) 
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TABLE I: Coefficients B , B 4 , B 3 , B 2 , B x for each model 



model 


Ni 


B 


B 4 


B 3 


B 2 


Bi 


Octagon 


8 


8 


48 


-6 


1 


l 

4 


Double Octagon 


16 


16 


24 


-3 


l 

2 


1 

8 


Triple Octagon 


24 


24 


16 


-2 


1 

3 


1 

12 



{£; + /} 



fc + / if k + / < 4 
fc + / - 4 if /c + / > 4 



And the coefficients £>o, 54 , i? 3 , 5 2 , Bi for each model are summarized in Table I. 

The temporal and spatial derivatives in Eq.(pQ) are computed by the FD schemes. Via 
the Chapman-Enskog expansion it is easy to find that the above FDLB model presents the 
same results as the following Navier-Stokes equations 

^ + V • (pu) = 0, (12) 
+ V • (PI + puu) + V • [MV • u)I - /x(V • u)I] = 0, (13) 

^(pE + X -pu 2 ) + V • [pu(E + l -u 2 + ?-)] 
-V • [k'VE + fxu ■ (Vu) - /m(V • u) + ^V« 2 ] = 0, (14) 

in the hydrodynamic limit, where ji and ft' are viscosity and heat conductivity, respectively. 
They are related to the pressure P and relaxation time r via 

k = 2/i = 2Pr. (15) 



B. Our contribution 



1. Operator- splitting scheme 

Equation (pD) could be written in the following scalar form in polar coordinates 

9fki . 9fki 1 dffci 1 r gg,.. 
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To present our new scheme, we first decompose the two-dimensional FDLB Eq. ( fT6l) into 
the following one-dimensional form 

dfki _ _1 [f _ feQ] 
dt ~ T [Jki Jkil 

^+Vkir?t=0 (17) 



dt 

9. 

dt 1 u to T dr 
dfki _l l v dfki n 



via the operator splitting scheme. 



2. Analytic solution for temporal evolution 

The first subequation in Eq. ( fTTI ) has a traditional discrete solution in the following form 

m At = fL-^(fL-f%) as) 

In fact the subequation can be given an analytical solution dynamically as below 

flt dt = exp(-|) [fL - m + m exp(|)] (19) 

Comparing with the Eq.ffT8l). Eq.ffT9l) gives an accurate solution, which is very helpful when 
the term ^ or (/^ — f^) is not small enough. 



3. MWB scheme for spatial evolution 



The last two subequations in Eq.(fT71) can be written uniformly as 



i> -> fki 

£ r or 9 
a -» v kir or -v kiQ . 



(20) 



Since v kir is a constant and r can also be regarded as a constant when consider the last 
subequation of Eq.ffTTl). we can further obtain 

^-a 2 ^-0 (21) 

By introducing the symbol, i/j(£j,t n ) = iff-, and performing the Taylor expansion, we get 

= J,? - aAt(^)] + ±a 2 At 2 (^)] + 0(At 3 ). (22) 
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The derivatives about £j in Eq. (l22l) are all calculated by using the second order upwind 
scheme, 



+ 0(A£ 2 ) if a>0 



_^!^ + 0(Aa if a<Q 



2A£ 



*?-W^-> +om tf a > 
v ,»_ 2 ^ +i+v ,» + 2 

A£ 2 



+ 0(Af) if a<0 
Thus, from Eq.(l22l we get the well-known Warming-Beam Scheme, 



(23) 



(24) 



n+1 



^ - - V>7_x) - |C(1 - C)(^ - 2^_x + V"- 2 ) if a > 
^ n - C(^ +1 - V;) + |C(1 + C)(^ - 2^ +1 + ^ +2 ) if a < 



(25) 



where the higher order tiny quantities have been omitted, C = aAt/Al; is the Courant- 
number. The Stability condition requires \C\ < 1. 

In this work we modify the Warming-Beam scheme. Firstly, Eq.f|25l) is written as 



n+1 _ „/,n 



1 



where 



^-[C + -C(1-\C\)(1-V)]S 



, ^ n -^_! if < C < 1 
' ^j+i-4>j if- 1 < c< 



(26) 



(27) 



if ^ ^ ^™ and a > 



const if = and a > 
^^g 1 if ^i^? anda<0 
const if = ^ and a < 
Then, by introducing a switch function S(r)), Eq.f|26l) is written as 



(28) 



^ n+1 = ^ - [C + 2^(1 - |C|)(1 - S{ V ))(1 ~ v)]6 
To make the scheme monotonous in space, we require 

< <p(C) < 1, 
where (p(C) is a quadratic polynomial function 



(29) 



(30) 



V(C) = \C\ + -\C\(1-\C\)(1-S( V ))(1- V ) 



(31) 



From Eqs.([30D-(|3U) we get 



p(0) = 
< <p(C) < 1 
V?(l) = p(-l) = 1. 



Equation ( 13T1) can be written as 



^(x) = x H — x(l — x)a 
2 

where a; = \C\, ip(x) = <p(C), a = (1 — — 77). Thus, Eq. flH^l) becomes 

^(0) = 
< ip{x) < 1 • 

V(i) = i. 

Equation (1331) describes a parabola which has a extremum value at 

1 1 

x e = - + -. 

2 a 

To satisfy all the three conditions in Eq.f|34l). we require 



x e > 1 or x e < 



From the conditions in Eq.f|36l) we have 



M < 2, 



i.e., 



So we choose 



\(l-S( V ))(l-r ] )\<2. 

sw - H - 1 



w + 1 

It should be pointed out that, besides the lattice Boltzmann equation, the 
works also for simulating hydrodynamic equations. 
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4- Combined scheme for the LB evolution 



By composing the solutions of the three subequations in Eq. (fT71h we get the combined 
scheme for the LB evolution, 



fir 1 = exp(-f )[f kl - m + r k ^Mf)} 

-[C r + \C r {l - \C r \)(l - S( Vr ))(l - Vr )]S r 
-[Co + \C 6 {1 - \C e \)(l - S( V o))(l - vo)}So 



where 



C r — Vki 



At 
* Ar 



C B 



r Ukie AO 



6f) = 



f ki (ir, i9) - f ki (ir - 1, i0) if v kir > 

f ki (ir + l,i0)- fki(ir : iO) if v kir < 

fki(ir, i0) - f k i(ir, i0 - 1) if v m > 

fki(ir, i0 + l)- fki(ir, i0) if v ki0 < 



Vr = < 



Ve = < 



fki(ir-l,i0)-f k i(ir-2,i0) 
fki(ir,i6)-f ki (ir-l,i6) 

const 

f M (ir+2,iO)-f M (ir+l,iO) 
f ki (ir+l,iO)-f ki (ir,iO) 

const 

fki{ir,iO-l)-f ki {ir,iQ-2) 
fki(ir,i0)-f ki (ir,i6-l) 

const 

fki{ir,iO+2)-f ki (ir,iQ+l) 
f ki (ir,iO+l)-f ki (ir,iO) 

const 



(39) 



(40) 



(41) 



(42) 



if fki(ir, iO) ^ f ki (ir - 1, iO) k v kir > 
if fki(ir, iO) = f ki {ir - 1, iO) k > 
if f ki (ir + 1, i0) ^ / fci (ir, i0) & w fcir < 
if /fci(«r + 1, i0) = / fci (ir, i0) k < 
if f ki {ir, iO) ^ f ki (ir, i9 - 1) k v kie > 
if f ki (ir, iO) = f ki (ir, i9 - 1) k v m > 
if f ki (ir, i9 + l)^ fki(ir, iO) k v kie < 
if fki(ir, iO) + 1 = / fci (ir, i0) k v U q < 
Here, ir and i0 are indexes of the coordinate. To this step, we have got a new conservative 
monotonous scheme with the second order accuracy. 



C. Boundary conditions 

The physical domain being considered in this work is an annular region whose radii are 
R2 > R\ > . When the inner radius R\ — > 0, the annular region approximates to a circular 
area. If the annular physical domain is periodic and the period is iVj along the circumferential 
direction, it can be divided into iVj parts, where iVj is just the total number of the directions 
of discrete velocity in the DVM. In this case, we can pick out only one part for calculations. 
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270 



FIG. 2: Sketch for the whole system and the computational domain with lattice nodes. 

If the period is N i /Nf ) we can pick out Nf connected parts as computational domain, where 
Nf is a positive integer. Thus, the computational domain is that with R\ < r < R 2 and 
< 9 < 2nNf/Ni. The computational domain has two boundaries in the radial direction and 
two in the circumferential direction. It is clear that periodic boundary conditions should be 
applied in the circumferential direction. On the other hand, the inner and outer boundaries 
in the radial direction should be treated specifically according to the specific situation under 
consideration. In this work we study the case with Ni — 8 and Nf — 1. Figure [2] shows a 
sketch for the computational domain. 

1. Radial boundary condition 

Assume that the total number of radial nodes is N r) radial increment is Ar = (i?2 — 
R\)/(N r — 1), and the radius r — R\ + (ir — 1) x Ar, where ir — 1, 2, • • • , N r . In the case 
where the densities can be considered continuous around the boundaries, we can obtain the 
density values on the ghost nodes (ir = —1, 0, iV r + 1, 7V r + 2) via linear interpolation scheme, 



or an interpolation scheme with the higher-order accuracy. The temperatures can be calcu- 
lated in a similar way. However, the determination of flow velocities depends on the specific 
situation under consideration. The simplest microscopic boundary condition is to assume 
that at each boundary node the system is in its thermodynamic equilibrium, i.e. fki = f^. 



p(ir, i6) = 2p(ir + 1, i 
p(ir, x9) = 2p(ir — 1, i 



9) - p(ir + 2,ii 
9) — p(ir — 2, i 



9) if ir < 1 

\0) if ir > N r 



(43) 
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(a) 



(b) 



FIG. 3: (a) Rotation from the azimuthal boundary with Ng to the one with 0. (b) Rotation from 
the azimuthal boundary with 1 to the one with Ng + 1. 

For the non-equilibrium microscopic boundary condition, the deviation from thermodynamic 
equilibrium, — f^, can be calculated via the interpolation scheme. 



2. Azimuthal boundary condition 

Similarly, assume that the total number of azimuthal nodes is Ng, azimuthal increment is 
A6 = 27rN f /(NiNg), and the angle 9 = id x A6, where id = 1, 2, • • • , N e . The distribution 
functions on ghost nodes (i9 — — 1, 0, Ng + 1, Ng + 2) are computed in the following way 



(44a) 



{i + Nf-Ni} 



f(ir,0,k,i) = f(ir,N e ,k,i + N f - N t ) 
f(ir, -1, k, i) = f(ir, Ng - 1, k,i + N f - N) 

i + Nf-Ni if i + N f > Ni 
i + Nf if i + Nf < N 

f(ir, N + 1, k, i) = f(ir, 1, k, i + N t - N f ) 
f(ir, Ng + 2, k, i) = f(ir, 2,k,i + N { - N f ) 

i-Nf + Ni if i - Nf < 
i - Nf if i - N f > ' 

A schematic diagram for the case with N — 8 and Nf = 1 is referred to FigJHl Figure (a) 
shows the way in which fki(ir, 0) is given from fki(ir, Ng) via rotation. Figure (b) shows the 
relation between fki(ir, Ng + 1) and fki(ir, 1) via rotation. 



(44b) 



{i-Nf + N} 
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III. VALIDATION AND VERIFICATION 



A. Performance on rotational slip flow 

We first consider the motion of a fluid between two coaxial cylinders, with radii Ri and 
i? 2 , rotating about their axis with angular velocities u\ and cl; 2 . It should be pointed out 
that, the compressibility of the fluid is proportional to the Mach number squared. Our 
PCLB model is for compressible fluid and physically consistent with this behavior. In this 
test the Mach number is small, so we roughly consider the system as incompressible. Due to 
the rotational symmetry, we have u r — 0, u$ = uq(t), P = P(r). For simplicity, we rewrite 
uq as u in this test. The Navier-Stokes equation for incompressible flow in cylindrical polar 
coordinates gives the following two equations: 

dP 

— pu 2 /r 0, (45a) 

or 

f d 2 u Idu u . ,aki \ 

+ = o. (45b) 

The second has the following solution, 

u = A + — (46) 
r 

where the constants A and B are found from the boundary conditions, 

R\uj\ for r — R\ 
R2UO2 for r = R 2 

As a result, we get the velocity distribution to be 

u 2 R 2 2 - ujiR\ R\RI(uji - oj 2 ) 1 

11-2 — ^ « 1 i<2 — ^1 

which is a non-slip Navier-Stokes solution and is independent of the viscosity \i. 

Initially, the system is in its thermodynamic equilibrium with p — 1.0, T — 1 and u r = 0, 
uq = 0. The other parameters are given as Ri — 1, R 2 — 2, uj\ = —0.1, W\ = 0.1, At = 10 -3 , 
A^ r x Nq — 100 x 20. Figure [4] shows the comparison of the LB results and the Analytic 
Solution(AS) for the final steady state, where the values of relaxation time r are 0.2, 0.1, 0.05, 
and 0.01, respectively. We can find that the simulation results have a satisfying agreement 
with AS. The slight mismatch is due to the weak compressibility of the fluid which is ignored 
in the analytical solution. 
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(47) 



0.20 




■ 


T=0.20 


• 


T=0.10 


A 


T=0.05 


T 


T=0.01 




AS 



1.4 1.6 
radius 



2.0 



FIG. 4: Comparison of LB results and analytical solution for the steady rotational velocity uq 
under various values of r. 




FIG. 5: The plane description of velocity field about the rotational slip flow at time t = 0.6: (a) the 
MWB scheme; (b) Lax-Wendroff scheme; (c) the second order up wind scheme; (d) Warming-Beam 
scheme. 

To compare our MBW scheme with other FD schemes for the convection term, we show, 
in Figj5j the simulation results of velocity field at time t = 0.6. The FD schemes used in 
Figs.(al)-(dl) are our MWB scheme, the Lax-Wendroff scheme, the second order upwind 
scheme, and the original Warming-Beam scheme, respectively. Figures (al)-(dl) show the 
global velocity fields in the computational domain. What Figs.(a2)-(d2) show are the en- 
largements of the portions in the corresponding squares in Figs.(al)-(dl). From Figj5](b2) 
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it is clear that the Lax-Wendroff scheme brings artificial oscillations in the tangential com- 
ponent of flow velocities near the inner and outer boundaries. Figures E](c2) and[5](d2) show 
that both the second order upwind and the original Warming-Beam schemes bring artificial 
oscillations in the radial component of flow velocities in the whole computational domain. 
The simulation result based on our new scheme has a satisfying agreement with theoretical 
analysis. 

B. Performance on discontinuity 

To check the performance of our new scheme on system with discontinuity, we consider a 
shock wave propagating outward in an annular system with radii Ri and R 2 . The physical 
quantities in shocked region and pre-shocked region satisfy the following Hugoniot relations, 

p(u - D) = p (u - D) 
< P-P = p {D - u ){u - u ) (49) 
E-E =I(P + P )(-L-I) 

where D is the velocity of shock wave. Due to the rotational symmetry, u means u r in the 
above equations. 

The initial physical field is below 

(p,Ur,u ,P) inner = (1.58824,0.785674,0,2.66667) , R r < r < R s 
{p,u r ,u Q ,P) outer = (1,0,0,1) , R s <r < R 2 

where R s is the position of shock front. We choose Ri = 2000, R 2 = 2025, Rs = 2001, 
D = 2.12132, r = 10- 5 , At = 10- 5 , N r x N e = 250 x 3. 

Figure [6] shows the simulation results of pressure P, density p, temperature T and velocity 
u r along the radius at time t — 8 using various schemes. From left to right, the four columns 
correspond to our MWB scheme, the Lax-Wendroff scheme, the second order upwind scheme, 
and the original Warming-Beam scheme, respectively. It is clear from the second column 
that the simulated physical fields from the Lax-Wendroff scheme have strong unphysical 
oscillations in the shocked region. The third column shows that the second order upwind 
scheme results in unphysical "overshoot" phenomena in the pressure, density, temperature 
and flow velocity at the shock front. The fourth column shows that the original Warming- 
Beam has the same drawback as that of the second order upwind scheme. Contrast to the 
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FIG. 6: The distributions of pressure P, density p, temperature T and velocity i£ r along the radius 
using various schemes. From left to right, the four columns correspond to the MWB scheme, 
Lax-Wendroff scheme, second order up wind scheme and Warming-Beam scheme. 

other three columns, the first column shows that the simulation results from our MWB are 
much more accurate and physically reasonable. 



C. Simulation study on Richtmyer-Meshkov instability 

The Richtmyer-Meshkov(RM) instability takes place when a shock wave travels across 
an interface separating two kinds of fluids. For the one-dimensional system with plane 
shock wave, several models have been proposed to describe the increase of the amplitude 
A. Roughly speaking, the increase rate of A first shows a linear relationship with itself, 
i.e., dA/dt = cA, where c is the increasing coefficient. In other words, the amplitude A 
increases exponentially with time according to the relation, A = A exp(ct). When the time 
t is very small, exp(ct) = 1 + ci, A = A + A ct. It is clear that, at the very beginning, 
the amplitude A linearly increases with time t. In the later time, the increasing coefficient 
c itself is no longer a constant any more. Therefore, the later stage is generally referred to 
as the nonlinear increasing stage. 

In 1960 Richtmyer (47j modified the linear theory of Taylor for Ray leigh- Taylor instability 
and proposed an impulsive model in the case of a reflected shock wave. The growth rate 
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reads, 

dA 7 . . . , , , Au. 

— = kAuA t A u A 1 = A (l - — ) 

where k = 2n/\ is the wave number, Au is the velocity change of the interface when shock 
passes, A t represents the post-shock Atwood number, A\ is the post-shock amplitude, A is 
i;he initial amplitude. Cmpr = 1 — Au/ D is defined as compression ratio. In 1969 Meshkov 
48| measured growth rate and found that it is only about one half of that predicted by the 
impulsive model. In 1992 Benjamin [49] got similar experiment results. In 1997 Zhang and 



Sohn 



50| developed the following model 

dA vq 



dt 1 + k 2 vM + max[0, (kA^ 2 - (A^ 2 + 0.5] (kv t) 2 

where v = kAuAiAi. This model is the growth rate of the interface amplitude from early 
to late times in the case of transition from light medium to heavy one. 

In an annular system with radii R\ — 1.0 and R2 = 2.0, we study the RM instability in 
the following two cases: shocking from light to heavy media and shocking from heavy to light 
media. We consider the case where an initial sinusoidal perturbation, r = R+A x cos(kR9) : 
is applied to the density field, where R is the average position of the interface between the 
two media. Even though the system considered here shows two-dimensional effects, for the 
case where the perturbation wave length is small and the inner radius is large enough, the 
above theory for one-dimensional system with plane shock wave still works approximately. 



1. Shocking from light to heavy media 

We consider the case where a shock wave travels outward from the light medium to the 
heavy one with the velocity D = 2. The initial physical field is given as below 

[p ) u r) ue ) P)inner = (1.5,0.666667,0,2.33333) 

< (p,U r ,Uo,P) middle = (1,1,0,1) 
(p,U r ,U ,P) outer = (3,1,0,1) 

V 

where the subscripts inner and middle indicate shocked and pre-shocked regions of light 
medium, outer represents outer region where the heavy medium locates. The numerical 
values between the shocked and pre-shocked regions satisfy with the Hugoniot relations. We 
choose A = 0.02, R = 1.2, k = 20, r = At = 10~ 5 , N r x N = 1000 x 450. 
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FIG. 7: Snapshots of RM instability for the case where the shock wave travels outwards from the 
light to heavy media. The density and pressure fields at the times, t = 0, t = 0.05, t = 0.5, t = 1.2, 
are shown from left to right, respectively. 

Figure [3 shows the snapshots of the density and pressure fields. The first row is for the 
density fields. The second is for the pressure fields. From left to right, the four columns are 
for the times, t — 0, t — 0.05, t = 0.5, t — 1.2, respectively. When the shock wave passes the 
interface, the perturbation amplitude A in the density field first decreases significantly due 
to compression of the shock wave, see Figs.(a-l)-(b-l), then it begins to increase under the 
pressure gradient, see Figs.(b-l)-(d-l), where asymmetric structures at the two sides of the 
interface eventually result in the occurrence of the bubbles in the light medium and spikes 
in the heavy medium. The corresponding pressure fields shown in Figs. (a-2)-(d-2) present 
complementary information for understanding the evolution of the density field. It should be 
pointed out that the misalignment of pressure and density gradients promotes deformation. 

In order to draw some quantitatively information to compare with the above theory, 
we show, from left to right in FigJEJ the amplitude, growth rate, average radial position 
and velocity of perturbation interface with time, where t = is defined as the time when 
the shock meets with the perturbation interface, and the amplitude is defined as half of 
the maximum distance from the crest to trough. From Fig. (a) we can observe clearly the 
compression and recovery of the perturbation amplitude. In the compressure process the 
amplitude drops rapidly, it decreases to the minimum A comp = 0.0143 at the time t comp = 0.03 
and recovers to its initial value at the time t reco = 0.13. For the convenience of analysis, 
a vertical line is plotted to guide the eyes in each of the Figs, (a)- (d) to indicate the time, 
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FIG. 8: Descriptions of RM instability for the case where shock wave travels outwards from light 
medium to heavy one. (a) Amplitude simulated by our model, (b) Growth rate of different models, 
(c) Averaged radial position, (d) Velocity of perturbed interface. A vertical line is shown in each 
plot to guide the eyes for the time when the interface width recovers to its original value. 

t reC) when the interface width recovers its original value. The evolution process of the 
perturbation amplitude before the time t rec is referred to as the initial compressure and 
recovery stage and after that is referred to as the increasing stage. From the minimum 
value in Fig. (a), we get the compression ratio Cmpr = 0.0143/0.02 = 0.715. However, the 
theoretical solution based on the initial conditions is Cmpr = 0.758. The deviation of the 
simulation result from the theoretical one is about 6%. In Fig.(b) the simulated growth rate 
of the perturbation amplitude A qualitatively agrees well with those theoretical models in the 
recover stage. It should be pointed out that, after taking into account the two-dimensional 
effects existing in this test system but are ignored by the theory, the simulation result shows 
a satisfying agreement with the theory. In Figs.(c) and (d) we compare the position and 
velocity of interface, respectively, for two cases. The first case is the simple one-dimensional 
problem where the plane shock wave passes the plane interface of two fluids. The second 
case is for the situation under consideration. We can find that the velocity of perturbation 
interface is slower than that for the case with plane shock wave. Physically, in the first 
case, the shock wave does not result in transverse flow velocity, the interface propagates 
in one direction at a constant velocity; while in the second case, two reasons result in the 
decreasing of the propagation velocity of the interface. On the one hand, the RM instability 
results in vortexes of flows. The transverse motion distracts part of the kinetic energy. 
Consequently, the propagation velocity in the original direction becomes slower. On the 
other hand, the geometrical structure of the setup leads to a slower propagation velocity. 
During the propagation process, the radius where the interface locates becomes larger, the 
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FIG. 9: Snapshots of RM instability for the case where the shock wave travels outwards from the 
heavy to light media. The density and pressure fields at the times, t = 0, 0.05, 0.5, and 1.2, are 
shown from left to right, respectively. 

area of the interface becomes larger. The distraction of kinetic energy makes slower the 
propagation velocity in the radial direction. 



2. Shocking from heavy to light media 

In the subsequent simulation, we choose Pi nner = 1, Pouter = 0.5, and other parameters 
are the same as those in the above simulation. 

We show the simulation results for density and pressure fields at times t = 0.0, 0.05, 
0.3, 0.5, and 1.2 in FigJHl where the interface reversal phenomenon is obviously observed. 
Simulation results show the following physical procedure: when the shock wave passes the 
interface, a rarefaction wave reflected inward and a transmission wave outward are generated. 
Near the crest, the pressure of heavy medium is smaller than that of the light medium. Driven 
by the pressure gradient, the perturbation amplitude decreases with the outward motion of 
the interface. Then, the peak and valley of initial interface invert, the heavy and light fluids 
gradually penetrate into each other as time goes on, the light fluid "falls" to form a bubble 
and the heavy fluid "rises" to generate a spike. 

We show the change of perturbation interface with time in FigHUl Figure (a) shows the 
simulated amplitude. Figure (b) shows the growth rate, where the line with scatters is for 
the LB result, the dashed line shows result from the Richtmyer model, and the dotted line is 
for the absolute value of result from the Richtmyer model. Figures (c)-(d) show the average 
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FIG. 10: Descriptions of RM instability for the case where shock wave travels outwards from 
heavy medium to light one. (a) Amplitude simulated by our model, (b) Growth rates from LB 
simulation, from the original Richtmyer model and from the absolute value of Richtmyer model, 
(c) Averaged radial position, (d) Velocity of perturbed interface. Three vertical lines are shown in 
each plot to guide the eyes for the compressure, reversal, recovery and increasing stages. 

radial position and velocity of perturbed interface. Three vertical lines are shown in Figs, (a) 
and (b) to divide the evolution progress into four stages, i.e., the stages of compression, 
reversal, recovery and increasing. The first guideline corresponds to the time, t comp = 0.03, 
when the amplitude is rapidly compressed to A comp = 0.013. The second one is for the time, 
tzero = 0.27, when the amplitude reaches zero. The last one indicates the time, t reco = 0.66, 
when the amplitude recovers to its initial value, A init = 0.02. From Fig. (a) we can get the 
compression ratio Cmpr = 0.013/0.02 = 0.65. As a comparison, the theoretical solution is 
Cmpr = 0.61. From Fig.(b) we can find that, for the evolution process of RM instability, 
our simulation result is close to that of Richtmyer model. Figures (c) and (d) show the same 
phenomena as those in FigJHl 



D. Simulation study on Kelvin-Helmholtz instability 

To investigate the Kelvin-Helmholtz (KH) instability in annular region with radii R\ < i?2, 
we set the initial physical field as below, 

/ \ Pinner H~ Pouter Pinner Pouter , T /T R\ /rn\ 

p(r) = tanh(-^— ), (50a) 

/ \ ^-inner ^-outer ^-inner ^-outer , ^ ^\ /rrvi \ 

u ( r ) = 2 2 tanh(-^— ), (50b) 

P(t) Pinner Pouter i (50c) 
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FIG. 11: Snapshots of KH instability for the case pinner < Pouter- The four columns are for the 
density and temperature contours at t = 0, 0.3, 0.7, and 1, respectively. 

where D p and D u are the widths of density and velocity transition layers. Pi nner , u inner and 
Pinner (Pouter, ^outer and Pouter) are the density, velocity and pressure of fluid near the inner 
(outer) cylinder, respectively. R is the radial position of interface between two media. In 
order to trigger the KH rollup, the following perturbation of velocity in the r-direction, 

u r e r = u e r sm(kR6) exp(— \r — (51) 

is added to the initial velocity field described by Eq. (l50bl) . where u = 0.5 is the amplitude 
of the perturbation, k is wave number. We study the KH instability in the following two 

Cases. Pinner ^ Pouter &Ild Pinner ^ Pouter- 



1. Case of Pinner ^ Pouter 

In the subsequent simulation, we choose p^ner = 0.5, p ou ter = 1-0, u outer = 0.5e^, u outer = 
-O.5e , D p = D u = 0.1, R x = 1, R 2 = 2, R = 1.5, k = 16, r = At = 10" 5 , N r x N = 
200 x 90. We plot the density and temperature contours at the times, t — 0, 0.3, 0.7, and 1, 
respectively in FigHH Figure (a) shows the initial state, where the interface starts to roll up 
gradually under the influence of the disturbance in Eq.(|5T1). In Figs.(b)-(d) we can observe 
that the deformation of interface caused by the KH instability become more significant with 
time. 

Let's study the spatial distribution of density and pressure at time t — 1 in Fig {121 Figures 
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FIG. 12: Snapshots at time £ = 1 for the case pinner < Pouter- ( a- l) an d (b-1) show the contours 
of density and pressure, (a-2) shows more clearly the contour of density with velocity field in the 
region labeled by the square in (a-l). (b-2) shows more clearly the contour of pressure with velocity 
field in the region labeled by the square in (b-1). 

(a-l) and (b-1) show the contours of density and pressure. Figure (a-2) shows more clearly 
the contour of density with velocity field in the region labeled by the square in Fig. (a-l). 
Figure (b-2) shows more clearly the contour of pressure with velocity field in the region 
labeled by the square in Fig. (b-1). From the velocity field in Fig. (a-2) we conceive that 
the KH instability would continue to develop and promote the intermixing and penetrating 
of the two fluids at the interface. It's clear to find in Fig. (b-2) that the pressure is in its 
minimum at the center of the vortex, and it is the pressure gradient that offers the centripetal 
force required for rotating. 
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FIG. 13: Snapshots of KH instability for the case pinner > Pouter- The four columns are for the 
density and temperature contours at t = 0, 0.3, 0.7, and 1, respectively. 

2. Case of p inner ^ Pouter 

In the subsequent simulation, pi nner = 1.0, p ou ter = 0.5, other parameters are the same 
as those in the case Pi nner < Pouter- We show the density and density contours at t — 0, 0.3, 
0.7, and 1, respectively in FigJT3l where the physical mechanism is similar to that in FigJTTl 
From Figs Hi] and [13] we find that the structures within the heavy medium are relatively 
sharp, likely "finger"; while the ones within the light medium are relatively smooth, likely 
"bubble". 

Figure [14] shows the spatial distributions of density, pressure and velocity at t — 1. From 
Figs JT2l and fT4l we find that the pressure is in its minimum at the vortex center and is in its 
maximum at the junction of vortices. 

It should be pointed out that, in the case of Pi nner > Pouter above, besides the KH 
instability, the Ray leigh- Taylor (RT) instability generally plays also a role in the rotating 
procedure of material. Because inertia presents an acceleration pointing to the light medium 
from the heavy one. But since the observation time is short, what we observe is mainly the 
result of the KH instability. 



25 




(a-2) (b-2) 

FIG. 14: Snapshots at time t = 1 for the case pinner > Pouter- ( a_ l) an d (b-1) show the contours 
of density and pressure, (a-2) shows more clearly the contour of density with velocity field in the 
region labeled by the square in (a-l). (b-2) shows more clearly the contour of pressure with velocity 
field in the region labeled by the square in (b-1). 

IV. NON-EQUILIBRIUM CHARACTERISTICS IN TWO SPECIFIC CASES 

To show the merit of the LB model over traditional ones, in this section we study the 
non-equilibrium characteristics in two specific cases. Among the seven moment relations, 
Eqs.([3l)-([9D, required by our model, only for the first three the equilibrium distribution 
function can be replaced by the distribution function f^. If we replace f%jt by fki in 
the left hand side of any one of Eqs.([6])-([9j), the left and right hand sides of Eqs.Q-Q will 
no longer be in balance. This mismatch measures the departure of the system from local 
thermodynamic equilibrium. 

We define two kinds of space-time dependent fields, moments M m and central moments 
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M^, as given below: 



< 



< 



M 2 (/fci) = ^2 ki fki^ki^ki 

M 3 (/^) = Ylki fki^ki^ki^ki 
M3,l(/fci) = Xlfci \fkiVki ' VkiVki 
< M4,2(/fci) = S/cz \fkiVki ' VkiVkiVki 

' M*(/ H ) = T,ki fki(v k i - u)(v fci - u) 

M 3(Ai) = E/cz Ai( v fci - u )( v /ci - U)(vjfei - U) 
M 3,l(/^) = Efci |/fc<(Vfei - U) • (v fci - U)(V W - U) 

k Ml >2 (f M ) = J2h \fki(y ki - u) • (v w - u)(v w - u)(v w - u) 



(53) 



(52) 



where the subscript "3, 1" means that the 3rd-order tensor is contracted to a lst-order 
tensor and the similar for "4,2". The moment M 35 i(= M^^^ a e a ) is a vector. It has two 
components, M^i^ and M^ifi. The moment M 2 (= M^^e^e^) is a second-order tensor 
with four components. Among the four components, only three, M^, M^rQ and M 2 ^ee ) are 
independent. The case for the moment 'WL^ 2 (= M^apGofip) is similar. The moment M 3 (= 
M 35a/ 3 7 e a e / 3e 7 ) is a third-order tensor with eight components. Among the eight components, 
only four, M 3?rrr , M 3?rr 0, M^^ee and M^^ee, are independent. The central moments is 
mathematically similar to M m . 

The central moment MJ; associates with the variance of the distribution function. Its trace 
associates with the temperature and its off-diagonal components associate with the shear 
effects. The former is a conserved quantity. When the system is not in its thermodynamic 
equilibrium state, the latter may not be zero. The moment M 3? i associates with the heat flux. 
For the equilibrium state, it describes convection of energy. For the non-equilibrium state, 
besides the energy convection, it includes also the thermodiffusion. For the central moment 
M3 l3 a breaking of the /(v) = /(— v) symmetry allows to eventually transport heat without 
necessarily carrying a net flow. The third-order central moment M3 may not be zero, while 
the first-order central moment MJ = ^2 ki /^(v^ — u) must be zero. In probability theory, for 
the one-dimensional distribution function /(v), the central moment M3 = f dvf(v)(v — u) 3 
is called "skewness". The fourth-order central moment Ml = J dvf(v)(v — u) 4 describes 
the "flatness" of the distribution and is also called "kurtosis". For a Gaussian distribution 
function, f(v) = 1/\/2tt exp[-(v - uf/2], Ml = 3. For the case with M 4 * > 3 and M 2 * = 1, 
the distribution is sharper than the Gaussian at the central position. 
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FIG. 15: Profiles of physical quantities (p, P, T, u r , uq) in the process of shock wave passing 
outwards from the heavy medium to the light one. The time t — 0.15. (a) Without initial 
perturbation at the interface, (b) With initial sinusoidal perturbation at the interface. Three lines 
are shown to guide the eyes for the interfaces. 

By Galilean invariance, it is clear that the moment M m contains the information of 
the macroscopic flow velocity u, while the moment MJ^ is only the manifestation of the 
thermo-fluctuations of molecules relative to the macroscopic flow velocity u. 

The manifestations of deviating from thermodynamic equilibrium from the two kinds of 
moments are as below: 

A m = M m (fki) ~ M m (/ fc e f) = M m (fki ~ ftd (54) 
K> = M* m (fki) ~ M* m (fJ) = M* m (fki ~ f%) (55) 

Similarly, A m contains the information of the macroscopic flow velocity u, while does 
not. 

A. Simulation results and analysis 

Now, we study the dynamic procedure where a shock wave propagates outwards from the 
heavy material to the light one. As the first step, we study the simplest situation where the 
incident shock wave is perpendicular to the unperturbed circular interface. In the second 
case, the interface is perturbed sinusoidally, and consequently the RM instability will occur. 
We choose such a time, t = 0.15, when the system shows three different interfaces. From 
left to right, the first is for the rarefaction wave, the second is for the two materials under 
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FIG. 16: Moments and their corresponding non-equilibrium manifestations for the case without 
perturbation at the material interface. The time t = 0.15. Figures (a)-(d) are for M2, M3, M3 5 i, 
M4 5 2, respectively. The symbols are for moments from and the solid lines are for moments 
from f^?. Figures (e)-(h) are for deviations A2, A3, A3 5 i, A^, respectively. Only independent 
components of M m and A m are shown. The specific correspondences are referred to the legends. 
Three squares are shown to guide the eyes for the interfaces. 

consideration, the third is for the shock wave. 

Figure [15] shows the profiles of physical quantities (p, P, T, u r) uq) along the radius 
with the azimuthal angle 6 = tt/6. Figure (a) corresponds to the case without initial 
perturbation at interface. Figure (b) is for the case with initial sinusoidal perturbation at 
the interface. Three lines are shown to guide the eyes for the interfaces. From Fig{T5] we 
can find that physical quantities change significantly at the three interfaces. For the case 
without perturbation in the material interface, we show the results of M m and A m in FigJT6l 
All independent components of M m , and A m are shown. The specific correspondences are 
referred to the legends. The 12 plots in Fig{T7| are the enlargements of the 12 portions 
labeled by the 12 squares in FigJTHl Figures HTlf a)-(d) correspond to the portions labeled by 
the first squares in Figs JTBlf e) -(h). respectively. Figures HTlf e)-(h) correspond to the portions 
labeled by the second squares in Figs JTBlf e) -(h). respectively. Figures HTlf i)-(l) correspond 
to the portions labeled by the third squares in Figs JTBlf e) -(h). respectively. The results of 
and are shown in FigJT8l The 12 plots in Fig{19] are the enlargements of the 12 
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FIG. 17: The deviations A m versus the radius, which are enlargements of the portions labeled by 
squares in Fig. [T6l Figures (a)-(d) are for the region around the first interface, with 1.00 < r < 1.15, 
in FigfT6l(e)-(h) respectively. Figures (e)-(h) are for the region around the second interface, with 
1.28 < r < 1.31. Figures (i)-(l) are for the region around the third interface, with 1.54 < r < 1.56. 

portions labeled by the 12 squares in FigUHl The specific correspondences between FigsJT9l 
and [18] are similar to the case of FigsH7] and [HI For the case with sinusoidal perturbation 
at the material interface, along the same radius, the results of M m and A m are shown in 
Figj20]and FigED The results of and A^ are shown in FigE2and Figj23l The specific 
correspondences between Figs ED and [20] and the specific correspondences between Figsl23l 
and f22l are also similar to the case of Figs JT71 and [T6l 

For both the two cases, one can clearly find the existence of the three interfaces via typical 
variations of the moments and corresponding moment differences. Around the interface for 
the shock wave, the system starts to deviate from thermodynamic equilibrium once the 
density starts to increase and goes back to its thermodynamic equilibrium as the density 
attains its steady value required by the Hugoniot relations. During the procedure deviating 
from thermodynamic equilibrium, A 2 , rr (or A^) shows a positive peak, while A 25 6><9 (or 
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FIG. 18: Central moments and their corresponding non-equilibrium manifestations for the case 
without perturbation at the material interface. The time t — 0.15. Figures (a)-(d) are for M|, M3, 
M3 l5 M4 2 1 respectively. The symbols are for central moments from and the solid lines are for 
central moments from f^. Figures (e)-(h) are for deviations A 2 , A3, A3 1? A \ 2 , respectively. Only 
independent components of and A^ are shown. The specific correspondences are referred to 
the legends. Three squares are shown to guide the eyes for the interfaces. 

A200) shows a negative peak with the same amplitude. For this interface, A 2j7 .0 or A 2r6 , is 
nearly zero. For the material interface, we can find two typical features. The peak value of 
A 2 ,rr or A 2 rr here is significantly smaller than that in the case with shock interface. The 
peak value of A 2 , r 6> or A 2r6 , is larger than that of A 25rr or A 2rr . The physical reason is as 
below. The shocking process is very quick, the shock interface is very thin. The changing 
rates of macroscopic quantities are very high. The system has very little time to relax 
to its thermodynamic equilibrium and has also very little time for the thermo diffusion. In 
contrast, due to the thermo diffusion with time, the material interface becomes wider and the 
system has enough time to relax to its thermodynamic equilibrium. During this procedure, 
the shear viscosity begins to play a role. Both the shock and rarefaction interfaces are 
mainly affected by mechanical effects, instead of the thermo diffusion. But the interface of 
rarefaction wave is much wider than that of the shock wave. Therefore, we can find that the 
system is much closer to its thermodynamic equilibrium around the rarefaction wave than 
around the shock wave. and A 4?2 (M4 2 and A£ 2 ) show similar behavior with M 2 and 
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FIG. 19: The deviations A^ versus the radius, which are enlargements of the portions labeled by 
squares in Fig. [TH1 Figures (a)-(d) are for the region around the first interface, with 1.00 < r < 1.15, 
in FigfT8l(e)-(h) respectively. Figures (e)-(h) are for the region around the second interface, with 
1.28 < r < 1.31. Figures (i)-(l) are for the region around the third interface, with 1.54 < r < 1.56. 



A 2 (Mg and Ag). Results of A3 and A3 x can be analyzed in a similar way. The components 
of A^ can be labeled by r p 6 q ) where p,q = 1,2, or 3. At the shock or rarefaction interface, 
if q — 0, the corresponding component is the largest. If q — 1 or 3, the corresponding 
component is negligibly small. 

Via comparing FigJTWa) and Figl2lTa) we can find that, the initial perturbation on the 
material interfaces enhances the shear viscosity effects. Other subfigures in FigJTTIand those 
in FigEU show consistent information. The information from A m in Fig{T7| (Figl2Tl) and 
that from A^ Fig{19] (Figl23l) are complimentary. 

Via comparing the interfaces for the shock and for the rarefaction in Figs JT6ll23l we can 
find that the shock wave increases density, pressure and temperature, while the rarefaction 
wave decreases those quantities. In other words, the two waves have opposite mechanical ef- 
fects. Although around both the two interfaces, from left to right, the density value becomes 
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FIG. 20: Moments and their corresponding non-equilibrium manifestations for the case with initial 
sinusoidal perturbation at the material interface. The time t — 0.15. Figures (a)-(d) are for M2, 
M3, M3 5 i, M4 5 2, respectively. The symbols are for moments from and the solid lines are for 
moments from f^. Figures (e)-(h) are for deviations A2, A3, A3 5 i, A^, respectively. Only 
independent components of M m and A m are shown. The specific correspondences are referred to 
the legends. Three squares are shown to guide the eyes for the interfaces. 



smaller, the non-equilibrium manifestations are oppositely different. The physical reason is 
as follows. The shock wave propagates outwards, while the rarefaction wave propagates 
inwards. Along their propagation directions, the temperature decreases around the shock 
wave, while it increases around the rarefaction wave. Via comparing the moments and cor- 
responding deviations for the two mechanical interfaces and those for the material interface, 
we can find that the latter shows a great difference in the non-equilibrium manifestations. 
The shear viscous effects are much more pronounced around the material interface. 

All the non-equilibrium effects in Figs JT6ll23l can be consistently interpreted as below. 
Among the four physical fields for the density, momentum, pressure and temperature, the 
temperature gradient is the most fundamental driving force triggering the non-equilibrium 
effects. The gradient of any other quantity triggers the non-equilibrium effects via triggering 
macroscopic transportation which leads to temperature gradient. The temperature gradient 
first initiates variance of the internal kinetic energy in the degree of freedom corresponding 
to the direction of the temperature gradient. (In this case the temperature shows gradient 
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FIG. 21: The deviations A m versus the radius, which are enlargements of the portions labeled by 
squares in Fig. [201 Figures (a)-(d) are for the region around the first interface, with 1.00 < r < 1.15, 
in Figl20l(e)-(h) respectively. Figures (e)-(h) are for the region around the second interface, with 
1.28 < r < 1.31. Figures (i)-(l) are for the region around the third interface, with 1.54 < r < 1.56. 



in the radial direction. This gradient first initiates the variance of the mean kinetic energy 
J dvf(v r — u r ) 2 /2.) Then, part of internal kinetic energy variance is transferred to other 
degrees of freedoms via collisions of molecules. Then, the internal kinetic energy in this 
degree of freedom further varies according to the temperature gradient, and so on. Only when 
the temperature gradient vanishes, the system can attain its thermodynamic equilibrium, 
i.e. the internal kinetic energy in different degrees of the freedom equal to each other. 



B. Recovering of the distribution function 

When the system is in its thermodynamic equilibrium state, the distribution function 
of the particle velocity is a local Maxellian, i.e. a normal distribution, symmetric about 
the mean flow velocity u, whose variance is proportional to the fluid temperature. These 
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FIG. 22: Central moments and their corresponding non-equilibrium manifestations for the case 
with initial sinusoidal perturbation at the material interface. The time t = 0.15. Figures (a)- 
(d) are for M^, M3, Mg l5 M42, respectively. The symbols are for central moments from 
and the solid lines are for central moments from f^. Figures (e)-(h) are for deviations A£, A3, 
Ag-L, A42, respectively. Only independent components of and A^ are shown. The specific 
correspondences are referred to the legends. Three squares are shown to guide the eyes for the 
interfaces. 

properties reflect profound symmetries of Newtonian mechanics, i.e. Galilean and scaling 
invariance, respectively. The local Maxwellian does not support any dissipative and trans- 
port mechanism, since these phenomena violate the aforementioned symmetries. Indeed, 
transport phenomena triggered by departures from local equilibria reflect into symmetry- 
breaking departures from the Maxwellian distribution. The maxwellian distribution is shown 
in FigJM 

From the results of the LB simulation, via the components of the deviations A^, we 
can draw qualitative information on the actual distribution function. As an example, we 
consider the case without perturbation at the material interface, and recover qualitatively 
the actual distribution function. The main steps are given below. 

We first consider the actual functions f(v r ) and f(vg) at the leftmost interface, the rar- 
efaction wave. It's easy to find in Fig JT9l( a) that A£ shows a negative peak and 
shows a positive peak with the same amplitude. Up to this step, we can imagine that the 
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FIG. 23: The deviations A^ versus the radius, which are enlargements of the portions labeled by 
squares in Fig. [22l Figures (a)-(d) are for the region around the first interface, with 1.00 < r < 1.15, 
in Figl22l(e)-(h) respectively. Figures (e)-(h) are for the region around the second interface, with 
1.28 < r < 1.31. Figures (i)-(l) are for the region around the third interface, with 1.54 < r < 1.56. 

distribution function f(v r ) is "thinner"and f(v$) is "fatter"than the Maxwellian. The peak 
of f(v r ) is higher and the peak of f(v$) is lower than that of the Maxwellian. A4 2 in FigJT9l 
(d) shows complementary information to A 2 in Fig JT9l( a). According to A3 in Fig JT9l( b) and 
A3 x in Fig JT9l( c). we can obtain that f(vo) is symmetric, while the f(v r ) is asymmetric. The 
portion for v r > is "fatter"than that for v r < 0. This is often called "positive skewness". A 
sketch of the actual distribution functions, f(v r ) and f{ve) ) are shown in Figj25] (a), where 
the dashed line shows the Maxwellian which is for both f eq (v r ) and f eq (vo). A sketch of the 
distribution functions around the shock wave is shown in Figl25l(c). where f(v r ) is "fatter"and 
f(vo) is a thinner"than the Maxwellian. The peak of f(v r ) is lower and the peak of f(vo) is 
higher than that of the Maxwellian and f(yo) is symmetric, while f(v r ) is asymmetric. The 
portion for v r > is "fatter"and the portion for v r < is "thinner". Similarly, we can obtain 
the sketch of actual distribution function around the materia interface. See Figl25l(b). 
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FIG. 24: The sketch of the Maxwellian distribution function in velocity space (v r ,vg). 



f(Vr) 



fist) 




(a) 



(b) 



(c) 



FIG. 25: The sketch of the Maxwellian and actual distribution functions versus velocity v r and 
vg, respectively. Figure (a) shows the distribution functions at the interface for rarefaction wave, 
Fig.(b) shows the ones at the material interface and Fig.(c) is for those at the interface for shock 
wave. The red solid line is for distribution function f(v r ), the blue solid one is for distribution 
function f(vg), and the dashed line is for f eq . 

Secondly, we study the contours of the actual distribution function in two-dimensional 
velocity space (v r ,vg). It's clear that the values of A?; rQ in Fig{T9] (a) and Fig{T9] (i) equal to 
zero, which implies that the contours of the actual distribution function at the rarefaction 
and the shock waves ought to be symmetric in v r or/and vg. With in mind that f(vg) is 
symmetric at the two interfaces, we can confirm that the v r coordinate axis is the symmetric 
axis for the contour. Figures [T9l(d) and (1) show consistent information. A^ r0 in FigJT9l 
(e) shows a positive peak, which implies that, at the material interface, the contour is not 
symmetric around the v r or vg coordinate. Because the shear viscous effects are pronounced, 
the actual distribution function is relatively complex. Figure l26l shows, from left to right, 
the sketches of contours of the actual distribution function at the interfaces of rarefaction, 
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(a) (b) (c) 

FIG. 26: The sketch of contours of the actual distribution functions in velocity space (v r ,v$). Figure 
(a)-(c) show the recovered distribution function contours at the three interfaces, rarefaction wave, 
material interface and shock wave, respectively. 




(a) (ft) ( c ) 



FIG. 27: The sketches of the actual distribution functions in velocity space (v r ,v$). Figures (a)- 
(c) show the recovered distribution functions at the three interfaces, rarefaction wave, material 
interface and shock wave, respectively. 

material and shock. 

Finally, by combining the results of the above two steps, we obtain the qualitative curves 
for the actual distribution functions at the three interfaces. The sketches are shown in 
Figl27l Figures (a)-(c) are for the rarefaction wave, the material interface and the shock 
wave, respectively. 
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V. CONCLUSIONS AND DISCUSSIONS 



A polar coordinate lattice Boltzmann model for compressible flows is presented. The key 
technique in this model is an operator-splitting scheme for the temporal and spatial evolu- 
tions. The convection term is solved via a modified Warming-Beam scheme where a switch 
function is introduced. The temporal evolution is calculated analytically. The new model 
works for both subsonic and supersonic flows, and consequently, it can be used to study com- 
plex flows under strong impact or shock. The new model is validated and verified via typical 
benchmark tests. For example, (i) the rotational slip flow, (ii) the Kelvin-Helmholtz(KH) 
instability, (iii) the stable shock tube problem, (iii) the Richtmyer-Meshkov(RM) instability. 
The latter two can not be simulated by the previous PCLB model j^j]. Even for the former 



two cases where the previous model 



451 ] works, the simulation results by the new model 



appear to be more accurate. Choosing computational domain and designing boundary con- 
ditions play an important role in numerical experiments. For annular systems showing 
periodic behaviors in the circumferential direction, one can pick out only one period of the 
domain for simulations. In such a case, the two boundaries in the circumferential direction 
are treated with periodic conditions. On the other hand, the two boundaries in the radial 
direction should be treated carefully according to the specific situation under investigation. 
The simplest microscopic radial boundary conditions assume that the system at the inner 
and outer boundaries are in thermodynamic equilibrium. More accurate distribution func- 
tions at the radial boundaries take into account the part deviating from equilibrium, which 
can be obtained via interpolation schemes from values at the neighboring lattice nodes inside 
the system. 

To show the merit of the PCLB model over the traditional methods based on hydrody- 
namics, we studied the macroscopic behaviors of the system resulted from departures from 
thermodynamic equilibrium around three kinds of interfaces, the shock wave, the rarefaction 
wave and the material interface, for two specific cases. In one of the two cases, the mate- 
rial interface is initially perturbed and consequently the RM instability occurs. It is found 
that, the macroscopic effects of deviating from thermodynamic equilibrium at the material 
interface are greatly different from those at the mechanical interfaces. The coupling effects 
of the molecular thermal motions in different degrees of freedom are much more pronounced 
at the material interface. The initial perturbation at the material interface results in more 
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pronounced two-dimensional effects and enhanced coupling effects of molecular motions in 
different degrees of freedom. The amplitude of deviating from equilibrium at the shock wave 
is much higher than those at the rarefaction wave and material interface. The above nu- 
merical results confirm that the temperature increase first initiates the increase of internal 
kinetic energy in the degree of freedom corresponding to the direction of temperature gra- 
dient, then increases those in other degrees of freedom, and vice verse. By comparing each 
component of the high-order moments and its value in equilibrium, we can draw qualitative 
information on the actual distribution function. These results deepen our understanding on 
the mechanical and material interfaces from a more fundamental point of view, and present 
valuable information to improve the physical modeling at the macroscopic level. Work based 
on the current LB model for systematic investigations of non-equilibrium effects in the RM 
instability and KH instability, is in progress. 

Acknowledgements 

The authors thank Prof. Guoxi Ni for many helpful discussions. AX and GZ acknowl- 
edge support of the Science Foundations of CAEP [under Grant Nos. 2012B0101014 and 
2011A0201002] and the Foundation of State Key Laboratory of Explosion Science and Tech- 
nology. AX, GZ, YL and CL acknowledge support of National Natural Science Foundation 
of China [under Grant Nos. 11075021, 11074300, and 91130020]. 



[1] S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University 

Press, New York, (2001). 

[2] F. J. Alexander, H. Chen, S. Chen and G. D. Doolen, Phys. Rev. A 46, 1967 (1992). 

[3] G. W. Yan, Y. S. Chen, S. X. Hu, Phys. Rev. E 59, 454 (1999). 

[4] C. H. Sun, Phys. Rev. E 58, 7283 (1998). 

[5] C. Sun and A. T. Hsu, Phys. Rev. E 68, 016303 (2003). 

[6] N. Cao, S. Chen, S. Jin, and D. Martinez, Phys. Rev. E 55, R21 (1997). 

[7] T. Kataoka and M. Tsutahara, Phys. Rev. E 69, 056702 (2004). 

[8] T. Kataoka and M. Tsutahara, Phys. Rev. E 69, 035701 (R) (2004). 



40 



[9] M. Watari and M. Tsutahara, Phys. Rev. E 67, 036306 (2003). 

[10] M. Watari and M. Tsutahara, Phys. Rev. E 70, 016703 (2004). 

[11] M. Watari, Physica A 382, 502 (2007). 

[12] A. Xu, Phys. Rev. E 71, 066706 (2005). 

[13] A. Xu, Europhys. Lett. 69, 214 (2005). 

[14] A. Xu, G. Gonnella, and A. Lamura, Phys. Rev. 67, 056105 (2003). 

[15] A. Xu, G. Gonnella, and A. Lamura, Phys. Rev. E 74, 011505 (2006). 

[16] A. Xu, G. Gonnella, A. Lamura, G. Amati, and F. Massaioli, Europhys. Lett. 71, 651 (2005). 

[17] A. Xu, G. Zhang, Y. Gan, F. Chen, and X. Yu, Front. Phys. 7(5), 582 (2012) 

[18] Y. Li, R. Shock, R. Zhang, and H. Chen, J. Fluid Mech. 519, 273 (2004). 

[19] X.F. Pan, A.G. Xu, G.C. Zhang, and S. Jiang, Int. J. Mod. Phys. C 18, 1747 (2007). 

[20] Y. Gan, A. Xu, G. Zhang, and Y. Li, Commun. Theor. Phys. 50, 201 (2008). 

[21] Y. Gan, A. Xu, G. Zhang, X. Yu, and Y. Li, Physica A 387, 1721 (2008). 

[22] Y. Gan, A. Xu, G. Zhang, and Y. Li, Commun. Theor. Phys. 56, 490 (2011). 

[23] Y. Gan, A. Xu, G. Zhang, and Y. Li, Phys. Rev. E 83, 056704 (2011). 

[24] Y. Gan, A. Xu, G. Zhang, Y. Li and H. Li, Phys. Rev. E 84, 046715 (2011). 

[25] Y. Gan, A. Xu, G. Zhang, and Y. Li, Europhys. Lett. 97, 44002 (2012). 

[26] Y. Gan, A. Xu, G. Zhang, and Y. Li, Front. Phys. 7(4), 481 (2012) 

[27] F. Chen, A. Xu, G. Zhang, Y. Gan, C. Tao, and Y. Li, Commun. Theor. Phys. 52, 681 (2009). 

[28] F. Chen, A. Xu, G.Zhang, Y. Li, S. Succi, EuroPhys. Lett. 90, 54003 (2010). 

[29] F. Chen, A. Xu, G. Zhang, Y. Li, Commun. Theor. Phys. 54, 1121, (2010). 

[30] F. Chen, A. Xu, G.Zhang, Y. Li, Commun. Theor. Phys. 55, 325 (2011). 

[31] F. Chen, A. Xu, G.Zhang, Y. Li, Phys. Lett. A 375, 2129 (2011). 

[32] F. Chen, A. Xu, G. Zhang, Y. Li, Commun. Theor. Phys. 56, 333, (2011). 

[33] F. Chen, A. Xu, G.Zhang, Y. Li, Theroe. & Appl. Mech. Lett. 1, 052004 (2011). 

[34] F. Nannelli, S. Succi, J. Stat. Phys. 68, 401 (1992). 

[35] Succi, Amati and Benzi, J. Stat. Phys. 81, 5 (1995). 

[36] Amati, Succi and Benzi, Fluid Dyn. Res. 19, 289 (1997). 

[37] G. Peng, H. Xi, C. Duncan and S. Chou, Phys. Rev. E 58, R4125 (1998). 

[38] G. Peng, H. Xi C. Duncan and S. Chou, Phys. Rev. E 59, 4676 (1999) 

[39] S. Ubertini, G. Bella and S. Succi, Phys. Rev. E 68, 016701 (2003) 

41 



[40] X. He, G. Doolen, J. Comput. Phys. 134, 306 (1997). 

[41] R. Mei, W. Shyy, J. Comput. Phys. 143, 426 (1998). 

[42] I. Halliday, L. A. Hammond, C. M. Care, K. Good, and A. Stevens, Phys. Rev. E 64, 011208 
(2001). 

[43] K. N. Premnath and J. Abraham, Phys. Rev. E 71, 056706 (2005). 

[44] P. Asinari, S. C. Mishra and R. Borchiellini, Numerical Heat Transfer B 57, 126 (2010). 

[45] M. Watari, Commun. Comput. Phys. 9, 1293 (2011). 

[46] P. Bhatnagar, E. P. Gross, and M. K. Krook, Phys. Rev. 94, 511 (1954). 

[47] R. D. Richtmyer, Comm. Pure Appl. Math. 13, 297 (1960). 

[48] E. E. Meshkov, Sov. Fluid Dyn. 4, 101 (1969). 

[49] R. F. Benjamin, Advances in Compressible Turbulent Mixing, edited by W. P. Dannevik, A. 

C. Buckingham, and C. E. Leith (1992). 

[50] Q. Zhang and S. Sohn, Phys. Fluids 9, 1106 (1997). 



42 



